
# -*- coding: mbcs -*-
import os
from part import *
from material import *
from section import *
from optimization import *
from assembly import *
from step import *
from interaction import *
from load import *
from mesh import *
from job import *
from sketch import *
from visualization import *
from connectorBehavior import *
# Modify the index to coordinate
os.chdir("./")
session.journalOptions.setValues(replayGeometry=COORDINATE, recoverGeometry=COORDINATE)
print("Modify the index to coordinate")

height,width,thickness=26,26,3
mdb.models['Model-1'].ConstrainedSketch(
    name='SKETCH',
    sheetSize=46)
mdb.models['Model-1'].sketches['SKETCH'].rectangle(
    point1=(0.0, 0.0), 
    point2=(height, width))
mdb.models['Model-1'].Part(
    dimensionality=THREE_D, 
    name='Part-1', 
    type=DEFORMABLE_BODY)
mdb.models['Model-1'].parts['Part-1'].BaseSolidExtrude(
    depth=thickness,
    sketch=mdb.models['Model-1'].sketches['SKETCH'])
del mdb.models['Model-1'].sketches['SKETCH']
print("Entity generate Successfully!")

mdb.models['Model-1'].Material(name='Material-1')
mdb.models['Model-1'].materials['Material-1'].Density(table=((4.56e-09, ), ))
mdb.models['Model-1'].materials['Material-1'].Elastic(table=((110000, 0.3), ))

mdb.models['Model-1'].materials['Material-1'].Plastic(
    hardening=JOHNSON_COOK, 
    scaleStress=None, 
    table=((1235, 1217, 0.2246, 1.3, 1080, 980), ))

# Section
mdb.models['Model-1'].HomogeneousSolidSection(
    material='Material-1', 
    name='Section-1', 
    thickness=None)
mdb.models['Model-1'].parts['Part-1'].SectionAssignment(
    offset=0.0,offsetField='', 
    offsetType=MIDDLE_SURFACE, 
    region=Region(
        cells=mdb.models['Model-1'].parts['Part-1'].cells.findAt(
            ((height/2.0, width/2.0, thickness), ), )), 
        sectionName='Section-1', 
        thicknessAssignment=FROM_SECTION)

mdb.models['Model-1'].rootAssembly.DatumCsysByDefault(CARTESIAN)
mdb.models['Model-1'].rootAssembly.Instance(
    dependent=ON, 
    name='Part-1-1', 
    part=mdb.models['Model-1'].parts['Part-1'])

mdb.models['Model-1'].ConstrainedSketch(
    gridSpacing=1.0, 
    name='SKETCH', 
    sheetSize=max(height,width)+40, 
    transform=mdb.models['Model-1'].parts['Part-1'].MakeSketchTransform(
        sketchPlane=mdb.models['Model-1'].parts['Part-1'].faces.findAt((height/2.0,width/2.0,thickness), ), 
        sketchPlaneSide=SIDE1,
        sketchUpEdge=mdb.models['Model-1'].parts['Part-1'].edges.findAt((height,width/2.0, thickness), ), 
        sketchOrientation=RIGHT, 
        origin=(height/2.0,width/2.0,thickness)))

#_infinite_position

infinite_thickness = 2
_height = height-2.0*infinite_thickness
_width = width-2.0*infinite_thickness

infinite_pt1 = (-height/2.0+infinite_thickness, width/2.0-infinite_thickness)
infinite_pt2 = (height/2.0-infinite_thickness, -width/2.0+infinite_thickness)

mdb.models['Model-1'].parts['Part-1'].projectReferencesOntoSketch(filter=
    COPLANAR_EDGES, sketch=mdb.models['Model-1'].sketches['SKETCH'])
mdb.models['Model-1'].sketches['SKETCH'].rectangle(point1=infinite_pt1, point2=infinite_pt2)
mdb.models['Model-1'].sketches['SKETCH'].Line(point1=(-height/2.0, width/2.0), point2=infinite_pt1)
mdb.models['Model-1'].sketches['SKETCH'].Line(point1=(_height/2.0, _width/2.0), point2=(height/2.0, width/2.0))
mdb.models['Model-1'].sketches['SKETCH'].Line(point1=(-_height/2.0, -_width/2.0), point2=(-height/2.0, -width/2.0))
mdb.models['Model-1'].sketches['SKETCH'].Line(point1=(_height/2.0, -_width/2.0), point2=(height/2.0, -width/2.0))
mdb.models['Model-1'].parts['Part-1'].PartitionFaceBySketch(
    faces=mdb.models['Model-1'].parts['Part-1'].faces.findAt(((height/2.0, width/2.0, thickness), )), 
    sketch=mdb.models['Model-1'].sketches['SKETCH'], 
    sketchUpEdge=mdb.models['Model-1'].parts['Part-1'].edges.findAt((height, width/2.0, thickness), ))
del mdb.models['Model-1'].sketches['SKETCH']

mdb.models['Model-1'].parts['Part-1'].PartitionCellByExtrudeEdge(
    cells=mdb.models['Model-1'].parts['Part-1'].cells.findAt(((height/2.0, width/2.0, thickness), )), 
    edges=(
        mdb.models['Model-1'].parts['Part-1'].edges.findAt((height/2.0, infinite_thickness, thickness), ),
        mdb.models['Model-1'].parts['Part-1'].edges.findAt((height/2.0, width-infinite_thickness, thickness), ), 
        mdb.models['Model-1'].parts['Part-1'].edges.findAt((infinite_thickness, width/2.0, thickness), ), 
        mdb.models['Model-1'].parts['Part-1'].edges.findAt((height-infinite_thickness, width/2.0, thickness), )), 
    line=mdb.models['Model-1'].parts['Part-1'].edges.findAt((0.0, 0.0, thickness/2.0), ), 
    sense=REVERSE)
mdb.models['Model-1'].parts['Part-1'].PartitionCellByExtrudeEdge(
    cells=mdb.models['Model-1'].parts['Part-1'].cells.findAt(((infinite_thickness/2.0, infinite_thickness/2.0, thickness), )), 
    edges=(
        mdb.models['Model-1'].parts['Part-1'].edges.findAt((height-infinite_thickness/2.0, infinite_thickness/2.0, thickness), ), 
        mdb.models['Model-1'].parts['Part-1'].edges.findAt((height-infinite_thickness/2.0, width-infinite_thickness/2.0, thickness), ), 
        mdb.models['Model-1'].parts['Part-1'].edges.findAt((infinite_thickness/2.0, infinite_thickness/2.0, thickness), ), 
        mdb.models['Model-1'].parts['Part-1'].edges.findAt((infinite_thickness/2.0, width-infinite_thickness/2.0, thickness), )), 
    line=mdb.models['Model-1'].parts['Part-1'].edges.findAt((0.0, 0.0, thickness/2.0), ), 
    sense=REVERSE)

datum_position = [[7.4, 7.4, 3], [7.4, 10.2, 3], [7.4, 13, 3], [7.4, 15.8, 3], [7.4, 18.6, 3], [10.2, 7.4, 3], [10.2, 10.2, 3], [10.2, 13, 3], [10.2, 15.8, 3], [10.2, 18.6, 3], [13, 7.4, 3], [13, 10.2, 3], [13, 13, 3], [13, 15.8, 3], [13, 18.6, 3], [15.8, 7.4, 3], [15.8, 10.2, 3], [15.8, 13, 3], [15.8, 15.8, 3], [15.8, 18.6, 3], [18.6, 7.4, 3], [18.6, 10.2, 3], [18.6, 13, 3], [18.6, 15.8, 3], [18.6, 18.6, 3], ]
mdb.models['Model-1'].rootAssembly.features.changeKey(
    fromName='Datum csys-1', 
    toName='Datum-csys-main')

for pos in datum_position:
    index = datum_position.index(pos)
    mdb.models['Model-1'].rootAssembly.DatumCsysByOffset(
        coordSysType=CARTESIAN, 
        datumCoordSys=mdb.models['Model-1'].rootAssembly.datums[1], 
        name='Datum csys-'+str(index+1), vector=(pos[0], pos[1], pos[2]))

for i in range(len(datum_position)):
    mdb.models['Model-1'].ExpressionField(
        description='', 
        expression='exp(-2*(sqrt(X**2+Y**2)/2.0)**10)', 
        localCsys=mdb.models['Model-1'].rootAssembly.datums[4+i], 
        name='AnalyticalField-'+str(i+1))

mdb.models['Model-1'].ExplicitDynamicsStep(
    improvedDtMethod=ON, 
    name='Step-1', 
    previous='Initial',
    timePeriod=3e-05)
for i in range(len(datum_position)-1):
    mdb.models['Model-1'].ExplicitDynamicsStep(
        improvedDtMethod=ON, 
        name='Step-'+str(i+2), 
        previous='Step-'+str(i+1),
        timePeriod=3e-05)
mdb.models['Model-1'].fieldOutputRequests['F-Output-1'].setValues(
    variables=("S","U","V","PE","PEEQ","ENER"))
# Amp
mdb.models['Model-1'].TabularAmplitude(
    data=((0, 0), (2e-09, 0.14), (4e-09, 0.26), (9e-09, 0.98), (1e-08, 1), (1.2e-08, 0.96), (1.8e-08, 0.8), (2.5e-08, 0.5), (8e-08, 0.2), (1.8e-07, 0.05), (2e-07, 0)),
    name='Amp-1', 
    smooth=SOLVER_DEFAULT, 
    timeSpan=STEP)

# Surf-Load
mdb.models['Model-1'].rootAssembly.Surface(
    name='Surf-Load', 
    side1Faces=mdb.models['Model-1'].rootAssembly.instances['Part-1-1'].faces.findAt(((width/2.0,height/2.0, thickness), )))
# Load
for i in range(len(datum_position)):
    mdb.models['Model-1'].Pressure(
        amplitude='Amp-1',

        createStepName='Step-'+str(i+1), 
        distributionType=FIELD, 
        field='AnalyticalField-'+str(i+1), 
        magnitude=3090.21, 
        name='Load-'+str(i+1), 
        region=mdb.models['Model-1'].rootAssembly.surfaces['Surf-Load'])
    if i<len(datum_position)-1: mdb.models['Model-1'].loads['Load-'+str(i+1)].deactivate('Step-'+str(i+2))    
# 2024-03-14 update


mdb.models['Model-1'].parts['Part-1'].setMeshControls(
    algorithm=ADVANCING_FRONT, 
    regions=mdb.models['Model-1'].parts['Part-1'].cells.findAt(
        ((infinite_thickness/2.0, width/2.0,thickness), ), 
        ((height/2.0, infinite_thickness/2.0, thickness), ), 
        ((height/2.0, width-infinite_thickness/2.0, thickness), ), 
        ((height-infinite_thickness/2.0, width/2.0, thickness), ), ), 
        technique=SWEEP)
mdb.models['Model-1'].parts['Part-1'].setSweepPath(
    edge=mdb.models['Model-1'].parts['Part-1'].edges.findAt((height-infinite_thickness/2.0, infinite_thickness/2.0, thickness), ), 
    region=mdb.models['Model-1'].parts['Part-1'].cells.findAt((height-infinite_thickness/2.0, width/2.0, thickness), ), 
    sense=FORWARD)
mdb.models['Model-1'].parts['Part-1'].setSweepPath(
    edge=mdb.models['Model-1'].parts['Part-1'].edges.findAt((height-infinite_thickness/2.0, infinite_thickness/2.0, thickness), ), 
    region=mdb.models['Model-1'].parts['Part-1'].cells.findAt((height/2.0, infinite_thickness/2.0,thickness), ), 
    sense=FORWARD)
mdb.models['Model-1'].parts['Part-1'].setSweepPath(
    edge=mdb.models['Model-1'].parts['Part-1'].edges.findAt((infinite_thickness/2.0, width-infinite_thickness/2.0, thickness), ), 
    region=mdb.models['Model-1'].parts['Part-1'].cells.findAt((height/2.0, width-infinite_thickness/2.0, thickness), ), 
    sense=REVERSE)
mdb.models['Model-1'].parts['Part-1'].setSweepPath(
    edge=mdb.models['Model-1'].parts['Part-1'].edges.findAt((infinite_thickness/2.0, width-infinite_thickness/2.0, thickness), ), 
    region=mdb.models['Model-1'].parts['Part-1'].cells.findAt((infinite_thickness/2.0, width/2.0, thickness), ), 
    sense=REVERSE)

# Mesh
mdb.models['Model-1'].parts['Part-1'].seedPart(
    deviationFactor=0.1,
    minSizeFactor=0.1, 
    size=0.15)

mdb.models['Model-1'].parts['Part-1'].seedEdgeByNumber(
    constraint=FINER,
    edges=mdb.models['Model-1'].parts['Part-1'].edges.findAt(
        ((infinite_thickness/2.0, infinite_thickness/2.0,thickness), ), 
        ((infinite_thickness/2.0, width-infinite_thickness/2.0, thickness), ), 
        ((height-infinite_thickness/2.0, infinite_thickness/2.0, thickness), ), 
        ((height-infinite_thickness/2.0, width-infinite_thickness/2.0, thickness), ), ), 
        number=1)
mdb.models['Model-1'].parts['Part-1'].seedEdgeByBias(
    biasMethod=SINGLE, 
    constraint=FINER, 
    end2Edges=mdb.models['Model-1'].parts['Part-1'].edges.findAt(
        ((0.0, 0.0, thickness/2.0), ), 
        ((height, 0.0, thickness/2.0), ), 
        ((height, width, thickness/2.0), ), 
        ((0.0, width, thickness/2.0), ), ), 
        maxSize=0.15, 
        minSize=0.05)
mdb.models['Model-1'].parts['Part-1'].generateMesh()

# Type
# C3D8R
mdb.models['Model-1'].parts['Part-1'].setElementType(
    elemTypes=(
        ElemType(
            elemCode=C3D8R, 
            elemLibrary=EXPLICIT, 
            secondOrderAccuracy=OFF, 
            kinematicSplit=AVERAGE_STRAIN, 
            hourglassControl=DEFAULT, 
            distortionControl=DEFAULT), 
        ElemType(
            elemCode=C3D6, 
            elemLibrary=EXPLICIT), 
        ElemType(
            elemCode=C3D4, 
            elemLibrary=EXPLICIT)), 
        regions=(mdb.models['Model-1'].parts['Part-1'].cells.findAt(
            ((height/2, width/2,thickness), )), 
        ))

mdb.models['Model-1'].parts['Part-1'].setElementType(
    elemTypes=(
        ElemType(
            elemCode=C3D8T, 
            elemLibrary=EXPLICIT, 
            secondOrderAccuracy=OFF, 
            distortionControl=DEFAULT), 
        ElemType(
            elemCode=C3D6T, 
            elemLibrary=EXPLICIT),
        ElemType(
            elemCode=C3D4T, 
            elemLibrary=EXPLICIT
        )), regions=(
        mdb.models['Model-1'].parts['Part-1'].cells.findAt(
            ((infinite_thickness/2.0, width/2.0,thickness), ), 
            ((height/2.0, infinite_thickness/2.0, thickness), ), 
            ((height/2.0, width-infinite_thickness/2.0, thickness), ), 
            ((height-infinite_thickness/2.0, width/2.0, thickness), ), ), ))

mdb.Job(
    activateLoadBalancing=False, 
    atTime=None, 
    contactPrint=OFF, 
    description='Do not submit the job', 
    echoPrint=OFF, 
    explicitPrecision=SINGLE, 
    historyPrint=OFF, 
    memory=90, 
    memoryUnits=PERCENTAGE, 
    model='Model-1', 
    modelPrint=OFF, 
    multiprocessingMode=DEFAULT, 
    name='DO_NOT_SUBMIT_THIS', 
    nodalOutputPrecision=SINGLE, 
    numCpus=2, 
    numDomains=2, 
    queue=None, 
    resultsFormat=ODB, 
    scratch='', 
    type=ANALYSIS, 
    userSubroutine='', 
    waitHours=0, 
    waitMinutes=0)

# Main Job
mdb.JobFromInputFile(
    activateLoadBalancing=False, 
    atTime=None, 
    explicitPrecision=SINGLE, 
    inputFileName='.//DO_NOT_SUBMIT_THIS.inp', 
    memory=90, 
    memoryUnits=PERCENTAGE, 
    multiprocessingMode=DEFAULT, 
    name="JOB_1_ANALYSIS", 
    nodalOutputPrecision=SINGLE, 
    numCpus=8, 
    numDomains=8, 
    queue=None, 
    resultsFormat=ODB, 
    scratch=r'./', 
    type=ANALYSIS, 
    userSubroutine='', 
    waitHours=0, 
    waitMinutes=0)

# Viewport
session.viewports['Viewport: 1'].setValues(displayedObject=mdb.models['Model-1'].parts['Part-1'])
session.viewports['Viewport: 1'].partDisplay.setValues(mesh=ON)
mdb.models['Model-1'].rootAssembly.regenerate()
print('''
             ******
             ******
             ******
             ******
             ******
             ******
             ******
             ******   
         ************
          ********** 
           ********
            ****** 
             ****  
              **
''')
print('**-Information-**: Please modify "C3D8T"->"CIN3D8" in "DO_NOT_SUBMIT_THIS.inp"file manually!')
